#include "./cell.hpp"
#include "yocto/code/utils.hpp"
#include "yocto/string/vfs-utils.hpp"

Cell:: ~Cell() throw()
{
}



////////////////////////////////////////////////////////////////////////////////
//
// Cell Initialization
//
////////////////////////////////////////////////////////////////////////////////

#define __CELL(field) field(Lua::Config::Get<lua_Number>(L,#field))

Cell:: Cell( const string &cfgile ) :
Chemistry(cfgile),
__CELL(surface),
__CELL(volume),
liters( 1e-15 * volume ),
__CELL(Cm),
Capa(Cm*surface*1e-14),
M( Library.size() ),
nvar( M + NEXTRA),
iEm( M + Em_OFFSET ),
sol( inside(0.0) ),
Em(0.0),
acc( Library ),
effectors(),
tmp( Library ),
perm_coeff( Library ),
ode_solver( 1e-4 ),
ode_derivs( this, & Cell::compute),
variables( nvar, 0.0 ),
__CELL(dt_max),
P_CO2( Lua::Function<double>(L,"P_CO2") )
{
    
    //==========================================================================
    // load the effectors
    //==========================================================================
    std::cerr << "-- Loading Effectors" << std::endl;
    const string cfgpath = _vfs::get_file_dir(cfgile);
    const string efffile = cfgpath + Lua::Config::Get<string>( L, "effectors" );
    Lua::Config::DoFile(L, efffile);
    _lua::load( L, effectors, "bio" );
    
    //==========================================================================
    // initialize solver
    //==========================================================================
    ode_solver.start(nvar);
    
    //==========================================================================
    // finalizing
    //==========================================================================
    std::cerr << "**** cell.surface         = " << surface << " microns^2" << std::endl;
    std::cerr << "**** cell.volume          = " << volume  << " microns^3" << std::endl;
    std::cerr << "**** surfacic capacity Cm = " << Cm << " microFarad/cm^2" << std::endl;
    std::cerr << "**** total capacity       = " << 1e12*Capa << " pF" << std::endl;
    for( size_t i=1 ; i <= M; ++i)
        perm_coeff(i) = 1.0;
}


void Cell:: ReloadInside()
{
    reload_inside(0.0);
    sol.copy( inside(0.0) );
}


void Cell:: load_state( const array<double> &data ) throw()
{
    assert( data.size() >= nvar);
    sol.get(data);
    Em = data[iEm];
}

void Cell:: save_state( array<double> &data ) const throw()
{
    assert( data.size() >= nvar);
    sol.put(data);
    data[iEm] = Em;
}


////////////////////////////////////////////////////////////////////////////////
//
// convert chemistry and effectors into an ODE
//
////////////////////////////////////////////////////////////////////////////////

void Cell:: compute(array<double> &dYdt, double t, const array<double> &Y)
{
    assert( nvar == dYdt.size() );
    assert( nvar == Y.size() );
    
    //--------------------------------------------------------------------------
    //-- load parameters
    //--------------------------------------------------------------------------
    const double    E   = Y[iEm];
    const solution &out = outside(t);
   
    tmp.get(Y); //-- temporary solution
    acc.ldz();  //-- reset accumulator
    
    //--------------------------------------------------------------------------
    //-- accumulator with leaks initialize leaks
    //--------------------------------------------------------------------------
    const double zeta = __Faraday__ * E / ( __R__ * Temperature );
    Leak( acc, zeta, tmp, out );
    
    //--------------------------------------------------------------------------
    //-- update accumulator with effectors
    //--------------------------------------------------------------------------
    for(Effectors::iterator i= effectors.begin(); i != effectors.end(); ++i )
    {
        const Effector &eff = **i;
        eff.call(L,acc,t,zeta,tmp);
    }
    
    
    //--------------------------------------------------------------------------
    //-- perform chemistry
    //--------------------------------------------------------------------------
    tmp.put( ChemSys.C  );
    acc.put( ChemSys.dC );
    ChemSys.reduce(t);
    acc.get( ChemSys.dC );
    
    
    //--------------------------------------------------------------------------
    // updated chemical part
    //--------------------------------------------------------------------------
    acc.put(dYdt);
    
    //--------------------------------------------------------------------------
    // update electrical part
    //--------------------------------------------------------------------------
    const double dQ = __Faraday__ * acc.sum_zC() * liters;
    //std::cerr << "dQ=" << dQ << " -> dE=" << dQ/Capa << std::endl;
    dYdt[iEm] = dQ / Capa;
}

void Cell:: Leak( solution &leak, double zeta, const solution &Sin, const solution &Sout ) const
{
    const double nrm = 1/volume;
    species::iterator         p  = Library.begin();
    component::iterator       i  = leak.begin();
    component::const_iterator j  = Sin.begin();
    component::const_iterator k  = Sout.begin();
    component::const_iterator q  = perm_coeff.begin();
    
    for( size_t count=M; count>0;--count,++p,++i,++j,++k,++q)
    {
        const species   &sp    = **p;
        Permeability    &Perm  = sp.get<Permeability>();
        component       &l     =  *i; assert( l.name    == sp.name );
        const component &Cin   =  *j; assert( Cin.name  == sp.name );
        const component &Cout  =  *k; assert( Cout.name == sp.name );
        const double     zz    = zeta * sp.z;
        const double     coeff = (*q).C;
        l.C += coeff * nrm * Perm.factor * Perm.SP(zeta) * Psi( zz ) * ( Cout.C - Cin.C * Exp( zz ) );
    }
    
}

void Cell:: Step( double t_ini, double t_end )
{
    //--------------------------------------------------------------------------
    // initialize variables
    //--------------------------------------------------------------------------
    save_state( variables );
    double  ode_h  = 1e-5;
    
    double t_curr = t_ini;
    while( t_curr < t_end )
    {
        const double t_next = min_of(t_curr + dt_max,t_end);
        
        //----------------------------------------------------------------------
        // perform step
        //----------------------------------------------------------------------
        ode_solver( ode_derivs, variables, t_curr, t_next, ode_h );
        
        //----------------------------------------------------------------------
        // forced normalization
        //----------------------------------------------------------------------
        t_curr = t_next;
        for( size_t i=M;i>0;--i)
            ChemSys.C[i] = variables[i];
        ChemSys.normalize(t_curr);
    }
    
    //--------------------------------------------------------------------------
    // return variables
    //--------------------------------------------------------------------------
    load_state( variables );
}


void Cell:: TryRescaleAllRates()
{
    lua_settop(L, 0);
    lua_getglobal(L, "rescale_rates" );
    if( lua_isnumber(L, -1) )
    {
        const double rescale_rates = lua_tonumber(L, -1);
        RescaleAllRates(rescale_rates);
    }    
}
